#coding:utf-8

from math import sqrt
import random
import numpy as np
import pylab as p
import mpl_toolkits.mplot3d.axes3d as p3

def gen_esferas(n):
	xi = [ 100*random.random() for j in range(0,n) ]
	yi = [ 100*random.random() for j in range(0,n) ]
	zi = [ 100*random.random() for j in range(0,n) ]
	
	ri = [ int( 100*random.random() ) for j in range(0,n) ]
	
	return (xi,yi,zi,ri)

n=3
res = gen_esferas(6)
xi = res[0]
yi = res[1]
zi = res[2]
ri = res[3]

'''
xi = [2,5,9]
yi = [2,5,9]
zi = [2,5,9]
ri = [5,7,1]

n = 3
'''

def soma_v(v1,v2):
	res = [0]*len(v1)
	for i in range(0,len(v1)):
		res[i] = v1[i]+v2[i]
	return res
	
def interno_v(v1,v2):
	res = 0
	for i in range(0,len(v1)):
		res += v1[i]*v2[i]
		
	return res	
	
def escalar_v(v1,k):
	res = [0]*len(v1)
	for i in range(0,len(v1)):
		res[i] = v1[i]*k
		
	return res	
	
def norma(v):
	res = 0
	for i in range(0,len(v)):
		res+=pow(v[i],2)
	res = sqrt(res)
	return res

	
def Ei(i,v):
	ind_i=3*i
	res = 0
	for k in range(0,n):
		ind_k=3*k			
		
		if i==k:
			continue
			
		parc = pow( (v[ind_i]-v[ind_k]),2 ) + pow( (v[ind_i+1]-v[ind_k+1]),2 ) + pow( (v[ind_i+2]-v[ind_k+2]),2 ) - pow( (ri[i]+ri[k]),2 )
		res+=pow(min(parc,0),2)
		
	return res
		
def teste(a,b,c,x,y,z,r1,r2):
	return pow((x-a),2)+pow((y-b),2)+pow((z-c),2)-pow((r1+r2),2)
	
def Fi(i,v):
	ind =3*i
	res = pow(v[ind],2)+pow(v[ind+1],2)+pow(v[ind+2],2) + Ei(i,v)
	#res = Ei(i,v)
	
	return res
	
def F(v):
	res = 0
	for i in range(0,n):
		res+=Fi(i,v)
		
	return res
	
def gen_dir(nvar):
	dirs = np.zeros( (nvar+1,nvar) )
	for i in range(0,nvar):
		dirs[i,i] = 1
		
	dirs[nvar,:] = [-1]*nvar
	
	return dirs
	
def main():
	delta = 2.0
	epson = 1e-4
	it = 0
	
	nvar = 3*n
	dirs = gen_dir(nvar)
	x0 = [0]*nvar
	
	f = F(x0)
	
	while (delta>epson and it<=6000):
		found_dir = False
		for d in dirs:
			xk = soma_v(x0,escalar_v(d,delta))
			f_ = F( xk )		
			
			if f_ < f:
				x0 = xk
				f=f_
				found_dir=True
				break
				
		if not(found_dir):
			delta/=2
			
		it+=1
	
	print it,delta
	return x0
		
def plot(x):
	u=np.r_[0:2*np.pi:100j]
	v=np.r_[0:2*np.pi:100j]

	px=np.ravel(np.cos(u)*np.sin(v))
	py=np.ravel(np.sin(u)*np.sin(v))
	pz=np.ravel(np.outer(np.ones(np.size(u)),np.cos(v)))
	
	for i in range(0,n):
		ind=3*i
		vx = np.array( [ x[ind] ]*len(px) )
		vy = np.array( [ x[ind+1] ]*len(py) )
		vz = np.array( [ x[ind+2] ]*len(pz) )
		
		(ri[i]*px+vx).tofile("x%d.txt" %(i),sep=' ')
		(ri[i]*py+vy).tofile("y%d.txt" %(i),sep=' ')
		z = (ri[i]*pz+vz)
		z.shape = (100,100)
		z.tofile("z%d.txt" %(i),sep=' ')
		
	
x = main()


#Print Scilab plot

print "v=[",
for i in range(0,len(x)):
	if i%3==0 and i!=0:
		print ";",
	print x[i],
print "]"

print "ri=",ri
